package com.ouc.eddyAlgorithm.utils;

import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

/**
 * 矢量几何方法探测涡旋（只适用北半球）
 */
public class VGEddyDetectUtil {

  private VGEddyDetectUtil() {
  }

  /**
   * 通过四个约束条件，返回区域涡旋矩阵
   * 返回一个二维矩阵。陆地（NaN）为-1，无涡旋为0，气旋为1，反气旋为2
   *
   * @param v           速度v分量
   * @param u           速度u分量
   * @param a           参数a
   * @param b           参数b
   * @param constraints 约束条件数组。数组中元素：true，表示添加该约束；false，表示不添加该约束
   * @return 符合约束条件的格网区域
   */
  public static int[][] getVgEddy(float[][] v, float[][] u, float[][] w, int a, int b, boolean[] constraints) {
    int m = v.length;
    int n = v[0].length;
    int[][] arr = new int[m][n];  //返回涡旋矩阵

    preEddyMatrixHandle(arr, w);  //涡旋矩阵预处理

    if (constraints[0]) { //第一约束条件
      firstConstraint(arr, v, a);
    }
    if (constraints[1]) { //第二约束条件
      secondConstraint(arr, u, a);
    }
    if (constraints[2]) { //第三约束条件
      thirdConstraint(arr, w, b);
    }
    if (constraints[3]) {  //第四约束条件
      fourConstraint(arr, u, v, a);
    }
    if (constraints[4]) {  //第四约束条件改进
      fourConstraintImp(arr, u, v, a);
    }

    return arr;
  }

  /**
   * 约束条件之前涡旋矩阵数据预处理
   * 对陆地进行预处理。陆地(NaN)设置为 -1，海洋设置为 0
   * 预处理步骤之后，矩阵数据只有-1和0
   *
   * @param arr 涡旋矩阵
   * @param w   流速
   */
  private static void preEddyMatrixHandle(int[][] arr, float[][] w) {
    int m = arr.length;
    int n = arr[0].length;

    //陆地设置为-1
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (Float.isNaN(w[i][j])) {
          arr[i][j] = -1;
        }
      }
    }

  }

  /**
   * 通过第一个约束条件，获取符合约束条件的区域范围
   * 修改涡旋矩阵。陆地（NaN）为-1，无涡旋为0，气旋为1，反气旋为2
   *
   * @param arr 涡旋矩阵
   * @param v   速度v分量
   * @param a   a 参数
   */
  private static void firstConstraint(int[][] arr, float[][] v, int a) {
    int m = v.length;
    int n = v[0].length;

    //设置气旋，反气旋
    for (int i = a + 1; i < m - a; i++) {
      for (int j = 1 + a; j < n - a; j++) {
        if (v[i][j] <= 0 && v[i][j - 1] >= 0) { //气旋（冷涡）
          if (v[i][j + a] < v[i][j] && v[i][j - a - 1] > v[i][j - 1]) {
            arr[i][j] = 1;
            arr[i][j - 1] = 1;
          }
        } else if (v[i][j] >= 0 && v[i][j - 1] <= 0) { // 反气旋（暖涡）
          if (v[i][j + a] > v[i][j] && v[i][j - a - 1] < v[i][j - 1]) {
            arr[i][j] = 2;
            arr[i][j - 1] = 2;
          }
        }
      }
    }
  }

  /**
   * 通过第二个约束条件，获取符合约束条件的区域范围
   * 修改涡旋矩阵。陆地（NaN）为-1，无涡旋为0，气旋为1，反气旋为2
   *
   * @param arr 经过第一个约束条件形成的涡旋矩阵区域
   * @param u   速度u分量
   * @param a   a 参数
   */
  private static void secondConstraint(int[][] arr, float[][] u, int a) {
    int m = arr.length;
    int n = arr[0].length;
    for (int i = a + 1; i < m - a; i++) {
      for (int j = 0; j < n; j++) {
        switch (arr[i][j]) {
          case 1: // 气旋（冷涡）
            if ((u[i][j] >= 0 && u[i - 1][j] <= 0) && (u[i + a][j] > u[i][j] && u[i - 1 - a][j] < u[i - 1][j])) {
              arr[i][j] = 1;
            } else {
              arr[i][j] = 0;
            }
            break;
          case 2: // 反气旋（暖涡）
            if ((u[i][j] <= 0 && u[i - 1][j] >= 0) && (u[i + a][j] < u[i][j] && u[i - 1 - a][j] > u[i - 1][j])) {
              arr[i][j] = 2;
            } else {
              arr[i][j] = 0;
            }
            break;
        }
      }
    }
  }

  /**
   * 通过第三个约束条件，获取符合约束条件的区域范围
   * 修改涡旋矩阵。陆地（NaN）为-1，无涡旋为0，气旋为1，反气旋为2
   *
   * @param arr 经过第二个约束条件形成的涡旋矩阵区域
   * @param w   流速
   * @param b   b 参数
   */
  private static void thirdConstraint(int[][] arr, float[][] w, int b) {
    int m = arr.length;
    int n = arr[0].length;
    for (int i = b; i < m - b; i++) {
      for (int j = b; j < n - b; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) {
          int regionCenterNum = arr[i][j];
          int[] regionMinIndex = findRegionMin(w, i, j, b);
          if (regionMinIndex[0] == -1) { //与局部最小值重合
            resetRegionMatrix(arr, i, j, regionCenterNum, b);
          } else {
            while (regionMinIndex[0] != -1) { //局部最小值不重合
              if (regionMinIndex[0] - b < 0 || regionMinIndex[0] + b >= m || regionMinIndex[1] - b < 0 || regionMinIndex[1] + b >= n) { //矩阵边界处理
                resetMatrixBoundaryNoEddy(arr, i, j, b);
                break;
              } else {
                resetRegionMatrix(arr, regionMinIndex[0], regionMinIndex[1], regionCenterNum, b);
                regionMinIndex = findRegionMin(w, regionMinIndex[0], regionMinIndex[1], b);
              }
            }
          }

        }
      }
    }

  }

  /**
   * 查找流速局部矩阵最小值
   * 局部矩阵若是最小值，返回[-1,-1];否则，返回最小值下标
   *
   * @param w 流速
   * @param i 局部矩阵中心
   * @param j 局部矩阵中心
   * @param b 参数 b
   * @return 返回最小值所在数组下标
   */
  private static int[] findRegionMin(float[][] w, int i, int j, int b) {
    float regionMin = w[i][j];
    int[] res = new int[]{-1, -1};
    for (int k = i - b; k <= i + b; k++) {
      for (int l = j - b; l <= j + b; l++) {
        if (w[k][l] < regionMin) {
          res[0] = k;
          res[1] = l;
        }
      }
    }
    return res;
  }

  /**
   * 以局部最小值为中心，重置局部矩阵
   * 设置局部矩阵中心为涡旋数据，其余为陆地或海洋
   *
   * @param arr             涡旋矩阵
   * @param i               局部矩阵中心
   * @param j               局部矩阵中心
   * @param regionCenterNum 局部矩阵中心数值
   * @param b               参数 b
   */
  private static void resetRegionMatrix(int[][] arr, int i, int j, int regionCenterNum, int b) {
    for (int k = i - b; k <= i + b; k++) {
      for (int l = j - b; l <= j + b; l++) {
        if (arr[k][l] == 1 || arr[k][l] == 2) {
          arr[k][l] = 0;
        }
      }
    }
    arr[i][j] = regionCenterNum;
  }

  /**
   * 重置边界局部矩阵。设为无涡旋局部矩阵
   * 应用于边界处理
   *
   * @param arr 涡旋矩阵
   * @param i   局部矩阵中心
   * @param j   局部矩阵中心
   * @param b   参数 b
   */
  private static void resetMatrixBoundaryNoEddy(int[][] arr, int i, int j, int b) {
    for (int k = Math.max(0, i - b); k < Math.min(arr.length, i + b); k++) {
      for (int l = Math.max(0, j - b); l < Math.min(arr[i].length, j + b); l++) {
        if (arr[k][l] == 1 || arr[k][l] == 2) {
          arr[k][l] = 0;
        }
      }
    }
  }

  /**
   * 通过第四个约束条件，获取符合约束条件的区域范围
   * 修改涡旋矩阵。陆地（NaN）为-1，无涡旋为0，气旋为1，反气旋为2
   *
   * @param arr 经过第三个约束条件形成的涡旋矩阵区域
   * @param u   水平速度
   * @param v   垂直速度
   * @param a   a 参数
   */
  private static void fourConstraint(int[][] arr, float[][] u, float[][] v, int a) {
    int m = arr.length;
    int n = arr[0].length;

    for (int i = a; i < m - a; i++) {
      for (int j = a; j < n - a; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) { //涡旋中心取值点
          if (!fourIsEddyCenter(u, v, a, i, j)) {
            arr[i][j] = 0;
          }
        }
      }
    }

  }

  /**
   * 应用于第四个约束条件，判断 arr[i][j] 点是否是涡旋中心.true:是涡旋中心；false：不是涡旋中心
   *
   * @param u 水平流速
   * @param v 垂直流速
   * @param a a 参数
   * @param i
   * @param j
   */
  private static boolean fourIsEddyCenter(float[][] u, float[][] v, int a, int i, int j) {
    int top = i - a + 1;  //上边界
    int bottom = i + a - 1; //下边界
    int left = j - a + 1; //左边界
    int right = j + a - 1;  //右边界
    //从下边界开始
    for (int k = left + 1; k < right; k++) {
      if (isOppositeBoundary(u[bottom][k - 1], v[bottom][k - 1], u[bottom][k], v[bottom][k])) {
        return false;
      }
    }
    //右边界
    for (int k = bottom - 1; k > top; k--) {
      if (isOppositeBoundary(u[k + 1][right], v[k + 1][right], u[k][right], v[k][right])) {
        return false;
      }
    }
    //上边界
    for (int k = right - 1; k > left; k--) {
      if (isOppositeBoundary(u[top][k + 1], v[top][k + 1], u[top][k], v[top][k])) {
        return false;
      }
    }
    //左边界
    for (int k = top + 1; k < bottom; k++) {
      if (isOppositeBoundary(u[k - 1][left], v[k - 1][left], u[k][left], v[k][left])) {
        return false;
      }
    }
    return true;
  }

  /**
   * 两点的流速方向是否是相反象限。true：对立象限；false：相同或相邻象限
   *
   * @param u1
   * @param v1
   * @param u2
   * @param v2
   * @return
   */
  private static boolean isOppositeBoundary(float u1, float v1, float u2, float v2) {
    if (u1 >= 0 && v1 >= 0) { //第一象限
      return u2 < 0 && v2 < 0;
    }
    if (u1 < 0 && v1 >= 0) { //第二象限
      return u2 >= 0 && v2 < 0;
    }
    if (u1 < 0 && v1 < 0) {  //第三象限
      return u2 >= 0 && v2 >= 0;
    }
    if (u1 >= 0 && v1 < 0) { //第四象限
      return u2 < 0 && v2 >= 0;
    }
    return false;
  }

  /**
   * 通过第四个约束条件改进，获取符合约束条件的区域范围
   * 改进之处在于：
   * 原第四约束条件：以涡旋中心为中心，长度为 2a-1 的矩阵边界进行方向判断，从而确定是否为涡旋中心；
   * 现第四约束条件改进：以涡旋中心为中心，长度从 1~2a-1 的多重矩阵进行方向判断，从而确定是否为涡旋矩阵。
   * 原第四约束条件可以解决的高涡度特征区域：1、流套区域；2、强反向流（相反流动的剪切区域）
   * 现第四约束条件改进可以进一步解决的问题：排除一定区域范围内多极值的情况。
   * 修改涡旋矩阵。陆地（NaN）为-1，无涡旋为0，气旋为1，反气旋为2
   *
   * @param arr 经过第三个约束条件形成的涡旋矩阵区域
   * @param u   水平速度
   * @param v   垂直速度
   * @param a   a 参数
   */
  private static void fourConstraintImp(int[][] arr, float[][] u, float[][] v, int a) {
    int m = arr.length;
    int n = arr[0].length;

    for (int i = a; i < m - a; i++) {
      for (int j = a; j < n - a; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) { //涡旋中心取值点
          if (!fourIsEddyCenterImp(u, v, a, i, j)) {
            arr[i][j] = 0;
          }
        }
      }
    }
  }

  /**
   * 应用于第四个约束条件改进，判断 arr[i][j] 点是否是涡旋中心.true:是涡旋中心；false：不是涡旋中心
   *
   * @param u 水平流速
   * @param v 垂直流速
   * @param a a 参数
   * @param i
   * @param j
   */
  private static boolean fourIsEddyCenterImp(float[][] u, float[][] v, int a, int i, int j) {
    //多重矩阵检测
    for (int k = 1; k <= a; k++) {
      if (!fourIsEddyCenter(u, v, k, i, j)) {
        return false;
      }
    }
    return true;
  }

  /**
   * 获取涡旋边界，边界是由八个方向点得出
   *
   * @param arr 涡旋中心点矩阵
   * @param v   速度v分量
   * @param u   速度u分量
   * @return 涡旋边界集合
   */
  public static List<int[]>[][] getVgEddyBoundary(int[][] arr, float[][] v, float[][] u) {
    int m = arr.length;
    int n = arr[0].length;
    List<int[]>[][] list = new ArrayList[m][n]; //返回涡旋边界集合

    int[] left, top, right, bottom;
    int[] leftTop, leftBottom, rightTop, rightBottom;  //八方向距离
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) {  //气旋
          list[i][j] = new ArrayList<>();
          left = getEddyBoundaryIndex(i, j, v, u, VgDirection.LEFT);
          top = getEddyBoundaryIndex(i, j, v, u, VgDirection.TOP);
          right = getEddyBoundaryIndex(i, j, v, u, VgDirection.RIGHT);
          bottom = getEddyBoundaryIndex(i, j, v, u, VgDirection.BOTTOM);
          leftTop = getEddyBoundaryIndex(i, j, v, u, VgDirection.LEFT_TOP);
          leftBottom = getEddyBoundaryIndex(i, j, v, u, VgDirection.LEFT_BOTTOM);
          rightTop = getEddyBoundaryIndex(i, j, v, u, VgDirection.RIGHT_TOP);
          rightBottom = getEddyBoundaryIndex(i, j, v, u, VgDirection.RIGHT_BOTTOM);
          Collections.addAll(list[i][j], left, top, right, bottom, leftTop, leftBottom, rightTop, rightBottom);
        }
      }
    }

    return list;
  }

  /**
   * 获取涡旋中心点半径
   * 沿东西方向，流速大小线性增加并方向相同，直至涡旋边界
   * 沿南北方向，流速大小线性增加并方向相同，直至涡旋边界
   * 获取四至方向至中心的半径大小，单位为格网点。
   * 取半径平均值，即为涡旋半径，保存至涡旋中心点位置
   *
   * @param arr 涡旋中心点矩阵
   * @param v   速度v分量
   * @param u   速度u分量
   * @return 涡旋中心点半径矩阵
   */
  public static float[][] getVgEddyRadius(int[][] arr, float[][] v, float[][] u) {
    int m = arr.length;
    int n = arr[0].length;
    float[][] radiusMatrix = new float[m][n]; //返回涡旋半径矩阵

    int left, top, right, bottom;
    int leftTop, leftBottom, rightTop, rightBottom;  //八方向距离
    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) {  //气旋
          left = j - getEddyBoundaryIndex(i, j, v, u, VgDirection.LEFT)[1];
          top = i - getEddyBoundaryIndex(i, j, v, u, VgDirection.TOP)[0];
          right = getEddyBoundaryIndex(i, j, v, u, VgDirection.RIGHT)[1] - j;
          bottom = getEddyBoundaryIndex(i, j, v, u, VgDirection.BOTTOM)[0] - i;
          int[] moveLeftTop = getEddyBoundaryIndex(i, j, v, u, VgDirection.LEFT_TOP);
          leftTop = (int) Math.sqrt(Math.pow(i - moveLeftTop[0], 2) + Math.pow(j - moveLeftTop[1], 2));
          int[] moveLeftBottom = getEddyBoundaryIndex(i, j, v, u, VgDirection.LEFT_BOTTOM);
          leftBottom = (int) Math.sqrt(Math.pow(i - moveLeftBottom[0], 2) + Math.pow(j - moveLeftBottom[1], 2));
          int[] moveRightTop = getEddyBoundaryIndex(i, j, v, u, VgDirection.RIGHT_TOP);
          rightTop = (int) Math.sqrt(Math.pow(i - moveRightTop[0], 2) + Math.pow(j - moveRightTop[1], 2));
          int[] moveRightBottom = getEddyBoundaryIndex(i, j, v, u, VgDirection.RIGHT_BOTTOM);
          rightBottom = (int) Math.sqrt(Math.pow(i - moveRightBottom[0], 2) + Math.pow(j - moveRightBottom[1], 2));
//          System.out.println(left + " " + top + " " + right + " " + bottom + " " + leftTop + " " + leftBottom + " " + rightTop + " " + rightBottom);
          radiusMatrix[i][j] = new BigDecimal(left + top + right + bottom + leftTop + leftBottom + rightTop + rightBottom).divide(new BigDecimal(8)).floatValue();
//          System.out.println("radiusMatrix[i][j] = " + radiusMatrix[i][j]);
        }
      }
    }

    return radiusMatrix;
  }

  /**
   * 获取涡旋中心点四至边界点坐标
   *
   * @param i           单个涡旋中心点坐标
   * @param j           单个涡旋中心点坐标
   * @param v           速度v分量
   * @param u           速度u分量
   * @param vgDirection
   * @return
   */
  private static int[] getEddyBoundaryIndex(int i, int j, float[][] v, float[][] u, VgDirection vgDirection) {
    int m = v.length;
    int n = v[0].length;

    switch (vgDirection) {
      case LEFT:
        while (j > 0 && Math.abs(v[i][j - 1]) > Math.abs(v[i][j])) {
          j--;
        }
        break;
      case TOP:
        while (i > 0 && Math.abs(u[i - 1][j]) > Math.abs(u[i][j])) {
          i--;
        }
        break;
      case RIGHT:
        while (j + 1 < n && Math.abs(v[i][j + 1]) > Math.abs(v[i][j])) {
          j++;
        }
        break;
      case BOTTOM:
        while (i + 1 < m && Math.abs(u[i + 1][j]) > Math.abs(u[i][j])) {
          i++;
        }
        break;
      case LEFT_TOP:
        while (i > 0 && j > 0 && Math.abs(v[i - 1][j - 1]) > Math.abs(v[i][j]) && Math.abs(u[i - 1][j - 1]) > Math.abs(u[i][j])) {
          i--;
          j--;
        }
        break;
      case LEFT_BOTTOM:
        while (i + 1 < m && j > 0 && Math.abs(v[i + 1][j - 1]) > Math.abs(v[i][j]) && Math.abs(u[i + 1][j - 1]) > Math.abs(u[i][j])) {
          i++;
          j--;
        }
        break;
      case RIGHT_TOP:
        while (i > 0 && j + 1 < n && Math.abs(v[i - 1][j + 1]) > Math.abs(v[i][j]) && Math.abs(u[i - 1][j + 1]) > Math.abs(u[i][j])) {
          i--;
          j++;
        }
        break;
      case RIGHT_BOTTOM:
        while (i + 1 < m && j + 1 < n && Math.abs(v[i + 1][j + 1]) > Math.abs(v[i][j]) && Math.abs(u[i + 1][j + 1]) > Math.abs(u[i][j])) {
          i++;
          j++;
        }
        break;
    }

    return new int[]{i, j};
  }

  /**
   * 获取涡旋中心点涡度ξ
   * ξ = dv/dx -du/dy
   *
   * @param arr 涡旋中心点矩阵
   * @param v   速度v分量
   * @param u   速度u分量
   * @param x   经度间隔
   * @param y   纬度间隔
   * @return 涡旋中心点涡度矩阵
   */
  public static float[][] getVgEddyAmp(int[][] arr, float[][] v, float[][] u, float x, float y) {
    int m = arr.length;
    int n = arr[0].length;
    float[][] ampMatrix = new float[m][n]; //返回涡旋涡度矩阵

    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) {
          ampMatrix[i][j] = (v[i][j] - v[i - 1][j]) / x - (u[i][j] - u[i - 1][j]) / y;
          ampMatrix[i][j] = Float.isNaN(ampMatrix[i][j]) ? 0.0f : ampMatrix[i][j];
        }
      }
    }

    return ampMatrix;
  }

  /**
   * 获取涡旋中心点涡动能（计算公式暂时存疑）
   * ξ = (u^2+v^2)/2
   *
   * @param arr 涡旋中心点矩阵
   * @param v   速度v分量
   * @param u   速度u分量
   * @return 涡旋中心点涡动能矩阵
   */
  public static float[][] getVgEddyEke(int[][] arr, float[][] v, float[][] u) {
    int m = arr.length;
    int n = arr[0].length;
    float[][] ekeMatrix = new float[m][n]; //返回涡旋涡动能矩阵
    int unit = 1000 * 100;  //转为cm

    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if (arr[i][j] == 1 || arr[i][j] == 2) {
          ekeMatrix[i][j] = (float) ((Math.pow(v[i][j], 2) + Math.pow(u[i][j], 2)) / 2 * unit);
        }
      }
    }

    return ekeMatrix;
  }

  /**
   * 获取涡旋中心点涡动能（计算公式暂时存疑）
   * ξ = (u^2+v^2)/2
   *
   * @param arr       涡旋中心点矩阵
   * @param ekeArr    涡动能
   * @param radiusArr 能量密度
   * @return 涡旋中心点涡动能矩阵
   */
  public static float[][] getVgEddyDens(int[][] arr, float[][] ekeArr, float[][] radiusArr) {
    int m = arr.length;
    int n = arr[0].length;
    float[][] densMatrix = new float[m][n]; //返回涡旋涡动能矩阵

    for (int i = 0; i < m; i++) {
      for (int j = 0; j < n; j++) {
        if ((arr[i][j] == 1 || arr[i][j] == 2) && radiusArr[i][j] != 0) {
          densMatrix[i][j] = BigDecimal.valueOf(ekeArr[i][j])
              .divide(BigDecimal.valueOf(Math.PI * Math.pow(radiusArr[i][j], 2)), 3, BigDecimal.ROUND_HALF_UP)
              .floatValue();
        }
      }
    }

    return densMatrix;
  }


}

enum VgDirection {
  LEFT, TOP, RIGHT, BOTTOM, LEFT_TOP, LEFT_BOTTOM, RIGHT_TOP, RIGHT_BOTTOM
}



